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Abstract 

We analyze the kinematics of multigrid Monte Carlo algorithms by inves- 
tigating acceptance rates for nonlocal Metropolis updates. With the help of a 
simple criterion we can decide whether or not a multigrid algorithm will have 
a chance to overcome critial slowing down for a given model. Our method is 
introduced in the context of spin models. A multigrid Monte Carlo procedure 
for nonabelian lattice gauge theory is described, and its kinematics is analyzed 
in detail. 
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15-19 September 1992, Amsterdam 



1 Introduction 

Multigrid Monte Carlo algorithms were introduced to overcome critical slowing down 
(CSD) []I|. For some interacting models, numerical experiments showed that the dy- 
namical critical exponent z was indeed substantially reduced p|. For other models, 
still z '^1 was found f^ . 

Therefore, an improved theoretical understanding of these algorithms is desir- 
able. We study the kinematics of multigrid Monte Carlo by investigating the scale 
dependence of acceptance rates for nonlocal Metropolis updates. A simple criterion 
for the possible success of a multigrid algorithm is given. Our analysis is also useful 
for the design of new multigrid algorithms for nonabelian gauge fields. 

2 Multigrid Monte Carlo 

First we consider models with partition functions 

Z= f Y[d<j)^ exp{-n{(j))) (1) 

on cubic d- dimensional lattices Aq. We shall use dimensionless spin variables. Non- 
local Monte Carlo updates are defined as follows: Divide the fundamental lattice Aq 
in cubic blocks of size l'^ (e.g. 1 = 2). This defines a block lattice Ai. By iterating 
this procedure one gets a hierarchy of block lattices Aq, Ai, . . . , Ak- We denote block 
lattice points in A^ by x'. Block spins ^x' are defined on block lattices A^. They 
are averages of the fundamental field over blocks of side length Lb = l^'- 

^x, = lV^"l-^''Y.<P^- (2) 

Xdx' 

A nonlocal change of the configuration (f) consists of a shift 

(px-^ <Px + SljJx, (3) 

where s is a real parameter. The shape of the nonlocal change is determined by the 
"coarse-to-fine interpolation kernel" if) that obeys the constraint 

Lb' E V'. = Lt'^^"^x',x'^ . (4) 

x(^x' 

Note that the effect of (|^) on A^ is $a./ — > ^x' + s for x' = x'^, whereas ^x' remains 
unchanged on the other blocks. The simplest t/^ is a piecewise constant kernel: 
ipx = Lb if x G x'g, and else. One can also use smooth kernels that avoid large 

energy costs from the block boundaries. 
The s-dependent Metropolis acceptance rate for such proposals is given by 

n{s) = (min[l, exp(-A7i:)]) , (5) 



where A7i = 7i(0 + sip) — H^cj)). The starting point of our acceptance analysis is 
the approximation formula [|, ^ ^^ 



Q{s) ^eiiciyhij . (6) 

Here, hi = (A7Y) denotes the average change in the fundamental Hamiltonian. 
Generally, the formula yields precise estimates that are confirmed by the acceptance 
rates directly measured in Monte Carlo simulations P, ^ . We will use (H) to predict 
the acceptance rate fl{s) for interacting models. Let us first discuss free massless 
field theory with action 7i(0) = ^{(j), — A0). Here, we obtain the exact result 



^(s)=erfc(^a/8|s|), (7) 

with a = (ipj—Atp). In d dimensions one finds a = 2dLB for piecewise constant 
kernels, and, for smooth kernels, a —* const if Lb ^^ 1. (For a systematic study 
of different kernels see ref. |@].) As a consequence, in massless free field theory, 
to maintain a constant acceptance rate (of, say, 50 percent) the amplitudes s have 

— 1/2 

to be scaled down like L^ for piecewise constant kernels, whereas for smooth 
kernels the acceptance rates do not depend on the block size. (At least for free field 
theory, the disadvantage of the piecewise constant kernels can be compensated for 
by using a W-cycle instead of a V-cycle. Smooth kernels can be used only in V-cycle 
algorithms.) 

3 Spin Models 

Now we discuss multigrid procedures for spin models. The kinematical analysis is a 

comparison of the scale dependence of acceptance rates for interacting models with 

the behavior in free field theory, where CSD is known to be eliminated by a multigrid 

algorithm. 

As an example we consider the 2-dimensional Sine Gordon model defined by the 

Hamiltonian 

^(0) = ^i^^-^^) - CEcos0.. (8) 

The model undergoes a (Kosterlitz-Thouless) phase transition at /3c, and Pc -^ Svr 
for ( ^ 0. In the massless phase (/? > (3c), the long wavelength excitations are spin 
waves. Since multigrid is an efficient method to accelerate spin waves in free massless 
field theory, one might expect that multigrid should be also the right method to fight 
CSD in the massless phase of the Sine Gordon model. But this expectation is wrong: 
For hi we find the expression 

hi = -s' + CCY.[^-cos{sij,)], (9) 

with C = (cos 0a;). Expanding for small s, we find that the second term in (P) 
behaves hke ~ s'^Y^^'^l: like a "mass" term in the kernel ip- Since J^xi^l is pro- 
portional to the block volume L^ for piecewise constant and for smooth kernels. 



one therefore has to face a dramatic decrease of acceptance when the blocks become 
large. A constant acceptance rate can only be achieved when the proposed steps 
s are scaled down like L~^^. It is therefore unlikely that any multigrid algorithm - 
based on nonlocal updates of the type discussed here - will be successful for this 
model. 

Analyzing multigrid algorithms, we found two classes of models P, |^. For 2d 
Sine Gordon and 0^ theory, s has to be rescaled like L^^ for piecewise constant 
and for smooth kernels, whereas for massless free field theory, the 2d XY model, 
the 2d 0{N) nonlinear a-model and U{1) lattice gauge theory, one can achieve 
L^-independent acceptance rates choosing smooth kernels. 

For almost all models of the second class, at least a substantial reduction of CSD 
could be achieved [0] . An exception is the 2d XY model in the vortex phase. There 
z ~ 1.4 was found 0. This shows that good acceptance rates alone are not sufficient 
to overcome CSD. 

The results of the analysis are consistent with the following rule: Sufficiently 
high acceptance rates for a complete elimination of CSD can only be expected if 
hi = (7i(0 + sip) — 7i(0)) contains no algorithmic "mass" term ~ s^ J2x 4^1- 

4 Nonabelian Gauge Fields 

We study 4-dimensional SU{2) lattice gauge theory with standard Wilson action 

mU) =/35:[l-iTrf/p] . (10) 

V 

The U-p are the usual path ordered products around plaquettes V. 

4.1 Covariant nonlocal update proposal 

We want to update only link variables Ux^r pointing in a selected direction r. Then 
only Ux,T variables in the same 3-dimensional slice A[ = {a; G Aq | x^ = t} are coupled 
(when all other link variables Ux,fj. with fi ^ r are kept fixed). Regard the Ux,t 
variables as spins. We then perform updates for a 3-dimensional spin model in a 
disordered background field given by the fixed link variables in the directions /x 7^ r. 
Now build 3-dimensional blocks x'^ of size L^ that are contained in the slice A[. 
Nonlocal update proposals will be global rotations of all link variables Ux,t attached 
to sites X inside the block x'^ by multiplying them with an S'f/(2)-matrix. 

Such global rotations only make sense if we have a certain smoothness of gauge 
fields inside the slice A[. We define a gauge transformation g by the Coulomb gauge 
condition 

Gc{U,g)= ^ Tr(^.[/,,^(?:+-) 1 max. (11) 

(x,x+A)eAj" 

If we actually performed this gauge transformation f/J'^ = gxUx^nQl+jx we would end 
up with the gauge fields fixed to the Coulomb gauge. However, we will use the 



transformation g only to achieve a smooth rotation of the Ux,t variables inside the 
block. This smooth rotation is defined as follows: 
Propose new link variables U'^ ^ by 

U,^r^K,r = 9:Rxg.U,,r, (12) 

with 

Rx{n, s) = cos{s4'x/2) + i sin(s'?/'^/2) n-a . (13) 

s is a uniformly distributed random number from the interval [—£,£:], ra is a 3- 
dimensional random unit vector, and a are the Pauli matrices, ip again denotes a 
kernel normalized as in (^). Finally calculate the associated change of the Hamilto- 
nian A?-^ and accept the proposed link variables with probability min[l, exp(— AH)]. 

One can show that this is a valid algorithm 0. For the detailled balance con- 
dition it is crucial to ensure that the updates are reversible. We have to get the 
same g before and after the move Ux^r -^ U'^^. This is fulfilled here, since only link 
variables Ux,^ with fi j^ t enter in the Coulomb gauge functional. 

Note that we do not have to extremize the gauge functional perfectly. If we 
always use the same algorithm (e.g. a given number of relaxation sweeps starting 
from g = 1), we will always get the same g and the nonlocal update is reversible. 



Table 1: Comparison of rriu with physical masses 



lattice size (3 mo \f^ w^o+ fno/ ^/^ ?^d/"^o+ 



16^ 2.4 0.4955(2) 0.258(2) 0.94(3) 1.92 0.53 

20^ 2.6 0.4650(2) 0.125(4) 0.52(3) 3.72 0.89 



Estimates for physical masses are taken from ref. 



4.2 Acceptance analysis for SU{2) 

Now we are going to check whether our proposal has a chance to overcome CSD. 
For the quantity hi = (A7^) that enters in the approximation formula (^ for the 
acceptance rates, we get (using piecewise constant ip) 

hi = 3A {Lb - l)Ll sm^sL^^^'^ /2) 

+ 6PL|[1-cos(sLb'/V2)], (14) 

with 

A = -^{i:r{{n-aUl^n-a -Ul^)Hi:^)) , (15) 

^x,^i = ^2:+A,T^x+f,/x^r,r5 ^^^^ P = /3/2(Trt/p). If we expand hi for small s, we see 
that the second term (coming from plaquettes that have one link in common with 
the block) is linear in Lb, the usual surface energy with piecewise constant kernels. 



But the first term (coming from plaquettes tliat are entirely inside tlie block) is 
again a "mass" term: it behaves like ~ As^L^. 

Let us discuss the origin of this term in more detail: In the weak coupling limit 
/3 — > cx) we have only pure gauges, i.e. U^^^ — >• 1 (with the gauge transformation g 
defined as above). Then the difference of the two terms contributing to A vanishes, 
and y4 — > 0. In this limit the unwanted mass term vanishes, and the acceptance 
rates behave just like in massless free field theory. 

For finite (3, however, [/| 7^ 1 and therefore A > because of the disorder of 
the gauge fields inside the block. Motivated by this, we identify the square root of 
A with a "disorder mass" m^i = yA. This disorder mass term will dominate the 
average energy change hi for large Lb 

h = (AH) ~ mls'Ll . (16) 

Close to criticality one might hope that rriD scaled with physical masses in the 
theory, say, the square root of the string tension n or the lowest glue ball mass mQ+ . 
Then, the fluctuations created by the algorithm would behave similar to the physical 
fluctuations, and we would expect the algorithm to behave well. But we have to 
face the fact that rriD will be dominated by the local disorder and decreases much 
slower than the physical masses in the system. This will be examined now. 

4.3 Monte Carlo study of mn 

We computed mD for several values of /3. To maximize Gc we used 50 Gauss-Seidel 
relaxation sweeps. Tests showed that increasing the number of relaxation sweeps 
beyond 50 hardly lowered mr) any further. 

In table 1 we display the ratios of the disorder mass rriD with physical masses. 
The results show that the disorder mass is nearly independent of (3 in the range 
studied, whereas the physical masses decrease by roughly a factor of two. In order 
to keep constant acceptance rate for increasing block size Lb, one again has to 
rescale s ~ I /Lb- Updating on large blocks gets essentially ineffective, and we have 
to expect CSD for this algorithm. 

4.4 Implementation and Test 

We implemented and tested the algorithm with 3-dimensional blocks using piecewise 
constant kernels and a W-cycle as described in ref. 0. We did not observe any 
substantial speed up compared to a local heat bath algorithm (not even a constant 
factor). This is in agreement with our prediction that updates on larger blocks are 
not efficient. 



Conclusions 

The kinematical mechanism that leads to a failure of multigrid algorithms is well 
described by our analysis. We hope that a better understanding of this problem can 
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lead to improved multigrid algorithms that can overcome kinematical obstructions 
stemming from an algorithmic "mass" term. 
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